home *** CD-ROM | disk | FTP | other *** search
/ AGA Toolkit '97 / The AGA Toolkit '97.iso / miscellaneous / science / maths / calc / source / zmath.h < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-07  |  14.1 KB  |  413 lines
  1. /*
  2.  * Copyright (c) 1994 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Data structure declarations for extended precision integer arithmetic.
  7.  * The assumption made is that a long is 32 bits and shorts are 16 bits,
  8.  * and longs must be addressible on word boundaries.
  9.  */
  10.  
  11. #ifndef    ZMATH_H
  12. #define    ZMATH_H
  13.  
  14. #include <stdio.h>
  15. #include "alloc.h"
  16. #include "endian.h"
  17. #include "longbits.h"
  18.  
  19. #include "have_stdlib.h"
  20. #ifdef HAVE_STDLIB_H
  21. # include <stdlib.h>
  22. #endif
  23.  
  24.  
  25. #ifndef ALLOCTEST
  26. # if defined(CALC_MALLOC)
  27. #  define freeh(p) (((void *)p == (void *)_zeroval_) ||            \
  28.             ((void *)p == (void *)_oneval_) || free((void *)p))
  29. # else
  30. #  define freeh(p) { if (((void *)p != (void *)_zeroval_) &&        \
  31.              ((void *)p != (void *)_oneval_)) free((void *)p); }
  32. # endif
  33. #endif
  34.  
  35.  
  36. typedef    int FLAG;            /* small value (e.g. comparison) */
  37. typedef int BOOL;            /* TRUE or FALSE value */
  38. typedef unsigned long HASH;        /* hash value */
  39.  
  40.  
  41. #if !defined(TRUE)
  42. #define    TRUE    ((BOOL) 1)            /* booleans */
  43. #endif
  44. #if !defined(FALSE)
  45. #define    FALSE    ((BOOL) 0)
  46. #endif
  47.  
  48.  
  49. /*
  50.  * NOTE: FULL must be twice the storage size of a HALF
  51.  *     LEN storage size must be <= FULL storage size
  52.  */
  53.  
  54. #if LONG_BITS == 64            /* for 64-bit machines */
  55. typedef unsigned int HALF;        /* unit of number storage */
  56. typedef int SHALF;            /* signed HALF */
  57. typedef unsigned long FULL;        /* double unit of number storage */
  58. typedef long LEN;            /* unit of length storage */
  59.  
  60. #define BASE    ((FULL) 4294967296)    /* base for calculations (2^32) */
  61. #define BASE1    ((FULL) (BASE - 1))    /* one less than base */
  62. #define BASEB    32            /* number of bits in base */
  63. #define    BASEDIG    10            /* number of digits in base */
  64. #define    MAXHALF    ((FULL) 0x7fffffff)    /* largest positive half value */
  65. #define    MAXFULL    ((FULL) 0x7fffffffffffffff) /* largest positive full value */
  66. #define MAXUFULL ((FULL)0xffffffffffffffff) /* largest unsigned full value */
  67. #define    TOPHALF    ((FULL) 0x80000000)    /* highest bit in half value */
  68. #define    TOPFULL    ((FULL) 0x8000000000000000)    /* highest bit in full value */
  69. #define MAXLEN    ((LEN)    0x7fffffffffffffff)    /* longest value allowed */
  70.  
  71. #else                    /* for 32-bit machines */
  72. typedef unsigned short HALF;        /* unit of number storage */
  73. typedef short SHALF;            /* signed HALF */
  74. typedef unsigned long FULL;        /* double unit of number storage */
  75. typedef long LEN;            /* unit of length storage */
  76.  
  77. #define BASE    ((FULL) 65536)        /* base for calculations (2^16) */
  78. #define BASE1    ((FULL) (BASE - 1))    /* one less than base */
  79. #define BASEB    16            /* number of bits in base */
  80. #define    BASEDIG    5            /* number of digits in base */
  81. #define    MAXHALF    ((FULL) 0x7fff)        /* largest positive half value */
  82. #define    MAXFULL    ((FULL) 0x7fffffff)    /* largest positive full value */
  83. #define MAXUFULL ((FULL)0xffffffff)    /* largest unsigned full value */
  84. #define    TOPHALF    ((FULL) 0x8000)        /* highest bit in half value */
  85. #define    TOPFULL    ((FULL) 0x80000000)    /* highest bit in full value */
  86. #define MAXLEN    ((LEN)    0x7fffffff)    /* longest value allowed */
  87. #endif
  88.  
  89. #define    MAXREDC    5            /* number of entries in REDC cache */
  90. #define    SQ_ALG2    20            /* size for alternative squaring */
  91. #define    MUL_ALG2 20            /* size for alternative multiply */
  92. #define    POW_ALG2 40            /* size for using REDC for powers */
  93. #define    REDC_ALG2 50            /* size for using alternative REDC */
  94.  
  95.  
  96. typedef union {
  97.     FULL    ivalue;
  98.     struct {
  99.         HALF Svalue1;
  100.         HALF Svalue2;
  101.     } sis;
  102. } SIUNION;
  103.  
  104.  
  105. #if !defined(BYTE_ORDER)
  106. #include <machine/endian.h>
  107. #endif
  108.  
  109. #if !defined(LITTLE_ENDIAN)
  110. #define LITTLE_ENDIAN    1234    /* Least Significant Byte first */
  111. #endif
  112. #if !defined(BIG_ENDIAN)
  113. #define BIG_ENDIAN    4321    /* Most Significant Byte first */
  114. #endif
  115. /* PDP_ENDIAN - LSB in word, MSW in long is not supported */
  116.  
  117. #if BYTE_ORDER == LITTLE_ENDIAN
  118. # define silow    sis.Svalue1    /* low order half of full value */
  119. # define sihigh    sis.Svalue2    /* high order half of full value */
  120. #else
  121. # if BYTE_ORDER == BIG_ENDIAN
  122. #  define silow    sis.Svalue2    /* low order half of full value */
  123. #  define sihigh sis.Svalue1    /* high order half of full value */
  124. # else
  125.    :@</*/>@:    BYTE_ORDER must be BIG_ENDIAN or LITTLE_ENDIAN    :@</*/>@:
  126. # endif
  127. #endif
  128.  
  129.  
  130. typedef struct {
  131.     HALF    *v;        /* pointer to array of values */
  132.     LEN    len;        /* number of values in array */
  133.     BOOL    sign;        /* sign, nonzero is negative */
  134. } ZVALUE;
  135.  
  136.  
  137.  
  138. /*
  139.  * Function prototypes for integer math routines.
  140.  */
  141. #if defined(__STDC__)
  142. #define MATH_PROTO(a) a
  143. #else
  144. #define MATH_PROTO(a) ()
  145. #endif
  146.  
  147. extern HALF * alloc MATH_PROTO((LEN len));
  148. #ifdef    ALLOCTEST
  149. extern void freeh MATH_PROTO((HALF *));
  150. #endif
  151.  
  152.  
  153. /*
  154.  * Input, output, and conversion routines.
  155.  */
  156. extern void zcopy MATH_PROTO((ZVALUE z, ZVALUE *res));
  157. extern void itoz MATH_PROTO((long i, ZVALUE *res));
  158. extern void utoz MATH_PROTO((FULL i, ZVALUE *res));
  159. extern void atoz MATH_PROTO((char *s, ZVALUE *res));
  160. extern long ztoi MATH_PROTO((ZVALUE z));
  161. extern FULL ztou MATH_PROTO((ZVALUE z));
  162. extern void zprintval MATH_PROTO((ZVALUE z, long decimals, long width));
  163. extern void zprintx MATH_PROTO((ZVALUE z, long width));
  164. extern void zprintb MATH_PROTO((ZVALUE z, long width));
  165. extern void zprinto MATH_PROTO((ZVALUE z, long width));
  166.  
  167.  
  168. /*
  169.  * Basic numeric routines.
  170.  */
  171. extern void zmuli MATH_PROTO((ZVALUE z, long n, ZVALUE *res));
  172. extern long zdivi MATH_PROTO((ZVALUE z, long n, ZVALUE *res));
  173. extern long zmodi MATH_PROTO((ZVALUE z, long n));
  174. extern void zadd MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  175. extern void zsub MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  176. extern void zmul MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  177. extern void zdiv MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res, ZVALUE *rem));
  178. extern void zquo MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  179. extern void zmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *rem));
  180. extern BOOL zdivides MATH_PROTO((ZVALUE z1, ZVALUE z2));
  181. extern void zor MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  182. extern void zand MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  183. extern void zxor MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  184. extern void zshift MATH_PROTO((ZVALUE z, long n, ZVALUE *res));
  185. extern void zsquare MATH_PROTO((ZVALUE z, ZVALUE *res));
  186. extern long zlowbit MATH_PROTO((ZVALUE z));
  187. extern long zhighbit MATH_PROTO((ZVALUE z));
  188. extern void zbitvalue MATH_PROTO((long n, ZVALUE *res));
  189. extern BOOL zisset MATH_PROTO((ZVALUE z, long n));
  190. extern BOOL zisonebit MATH_PROTO((ZVALUE z));
  191. extern BOOL zisallbits MATH_PROTO((ZVALUE z));
  192. extern FLAG ztest MATH_PROTO((ZVALUE z));
  193. extern FLAG zrel MATH_PROTO((ZVALUE z1, ZVALUE z2));
  194. extern BOOL zcmp MATH_PROTO((ZVALUE z1, ZVALUE z2));
  195.  
  196.  
  197. /*
  198.  * More complicated numeric functions.
  199.  */
  200. extern long iigcd MATH_PROTO((long i1, long i2));
  201. extern void zgcd MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  202. extern void zlcm MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  203. extern void zreduce MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *z1res, ZVALUE *z2res));
  204. extern void zfact MATH_PROTO((ZVALUE z, ZVALUE *dest));
  205. extern void zlcmfact MATH_PROTO((ZVALUE z, ZVALUE *dest));
  206. extern void zperm MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  207. extern void zcomb MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  208. extern FLAG zjacobi MATH_PROTO((ZVALUE z1, ZVALUE z2));
  209. extern void zfib MATH_PROTO((ZVALUE z, ZVALUE *res));
  210. extern void zpowi MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  211. extern void ztenpow MATH_PROTO((long power, ZVALUE *res));
  212. extern void zpowermod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
  213. extern BOOL zmodinv MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  214. extern BOOL zrelprime MATH_PROTO((ZVALUE z1, ZVALUE z2));
  215. extern long zlog MATH_PROTO((ZVALUE z1, ZVALUE z2));
  216. extern long zlog10 MATH_PROTO((ZVALUE z));
  217. extern long zdivcount MATH_PROTO((ZVALUE z1, ZVALUE z2));
  218. extern long zfacrem MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *rem));
  219. extern void zgcdrem MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  220. extern long zdigits MATH_PROTO((ZVALUE z1));
  221. extern FLAG zdigit MATH_PROTO((ZVALUE z1, long n));
  222. extern BOOL zsqrt MATH_PROTO((ZVALUE z1, ZVALUE *dest));
  223. extern void zroot MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *dest));
  224. extern BOOL zissquare MATH_PROTO((ZVALUE z));
  225. extern HASH zhash MATH_PROTO((ZVALUE z));
  226.  
  227.  
  228. /*
  229.  * Prime related functions.
  230.  */
  231. extern void zpfact MATH_PROTO((ZVALUE z, ZVALUE *dest));
  232. extern BOOL zprimetest MATH_PROTO((ZVALUE z, long count));
  233. extern long zlowfactor MATH_PROTO((ZVALUE z, long count));
  234.  
  235.  
  236. /*
  237.  * Misc misc functions. :-)
  238.  */
  239. #if 0
  240. extern void zapprox MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res1, ZVALUE *res2));
  241. extern void zmulmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
  242. extern void zsquaremod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  243. extern void zsubmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
  244. #endif
  245. extern void zminmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
  246. extern BOOL zcmpmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3));
  247.  
  248.  
  249. /*
  250.  * These functions are for internal use only.
  251.  */
  252. extern void ztrim MATH_PROTO((ZVALUE *z));
  253. extern void zshiftr MATH_PROTO((ZVALUE z, long n));
  254. extern void zshiftl MATH_PROTO((ZVALUE z, long n));
  255. extern HALF *zalloctemp MATH_PROTO((LEN len));
  256. extern void initmasks MATH_PROTO((void));
  257.  
  258.  
  259. /*
  260.  * Modulo arithmetic definitions.
  261.  * Structure holding state of REDC initialization.
  262.  * Multiple instances of this structure can be used allowing
  263.  * calculations with more than one modulus at the same time.
  264.  * Len of zero means the structure is not initialized.
  265.  */
  266. typedef    struct {
  267.     LEN len;        /* number of words in binary modulus */
  268.     ZVALUE mod;        /* modulus REDC is computing with */
  269.     ZVALUE inv;        /* inverse of modulus in binary modulus */
  270.     ZVALUE one;        /* REDC format for the number 1 */
  271. } REDC;
  272.  
  273. extern REDC *zredcalloc MATH_PROTO((ZVALUE z1));
  274. extern void zredcfree MATH_PROTO((REDC *rp));
  275. extern void zredcencode MATH_PROTO((REDC *rp, ZVALUE z1, ZVALUE *res));
  276. extern void zredcdecode MATH_PROTO((REDC *rp, ZVALUE z1, ZVALUE *res));
  277. extern void zredcmul MATH_PROTO((REDC *rp, ZVALUE z1, ZVALUE z2, ZVALUE *res));
  278. extern void zredcsquare MATH_PROTO((REDC *rp, ZVALUE z1, ZVALUE *res));
  279. extern void zredcpower MATH_PROTO((REDC *rp, ZVALUE z1, ZVALUE z2, ZVALUE *res));
  280.  
  281.  
  282. /*
  283.  * macro expansions to speed this thing up
  284.  */
  285. #define ziseven(z)    (!(*(z).v & 01))
  286. #define zisodd(z)    (*(z).v & 01)
  287. #define ziszero(z)    ((*(z).v == 0) && ((z).len == 1))
  288. #define zisneg(z)    ((z).sign)
  289. #define zispos(z)    (((z).sign == 0) && (*(z).v || ((z).len > 1)))
  290. #define zisunit(z)    ((*(z).v == 1) && ((z).len == 1))
  291. #define zisone(z)    ((*(z).v == 1) && ((z).len == 1) && !(z).sign)
  292. #define zisnegone(z)    ((*(z).v == 1) && ((z).len == 1) && (z).sign)
  293. #define zistwo(z)    ((*(z).v == 2) && ((z).len == 1) && !(z).sign)
  294. #define zisleone(z)    ((*(z).v <= 1) && ((z).len == 1))
  295. #define zistiny(z)    ((z).len == 1)
  296.  
  297. /*
  298.  * zgtmaxfull(z)    TRUE if abs(z) > MAXFULL
  299.  */
  300. #define zgtmaxfull(z)    (((z).len > 2) || (((z).len == 2) && (((SHALF)(z).v[1]) < 0)))
  301.  
  302. /*
  303.  * Some algorithms testing for values of a certain length.  Macros such as
  304.  * zistiny() do this well.  In other cases algorthms require tests for values
  305.  * in comparison to a given power of 2.  In the later case, zistiny() compares
  306.  * against a different power of 2 on a 64 bit machine.  The macros below
  307.  * provide a tests against powers of 2 that are independent of the work size.
  308.  *
  309.  *     zge16b(z)    TRUE if abs(z) >= 2^16
  310.  *     zge24b(z)    TRUE if abs(z) >= 2^24
  311.  *     zge31b(z)    TRUE if abs(z) >= 2^31
  312.  *     zge32b(z)    TRUE if abs(z) >= 2^32
  313.  */
  314. #if LONG_BITS == 64            /* for 64-bit machines */
  315. #define zge16b(z)    (!zistiny(z) && ((z).v[0] >= (HALF)0x10000))
  316. #define zge24b(z)    (!zistiny(z) && ((z).v[0] >= (HALF)0x1000000))
  317. #define zge31b(z)    (!zistiny(z) || (((SHALF)(z).v[0]) < 0))
  318. #define zge32b(z)    (!zistiny(z))
  319. #else
  320. #define zge16b(z)    (!zistiny(z))
  321. #define zge24b(z)    (((z).len > 2) || (((z).len == 2) && ((z).v[1] >= (HALF)0x100)))
  322. #define zge31b(z)    (((z).len > 2) || (((z).len == 2) && (((SHALF)(z).v[1]) < 0)))
  323. #define zge32b(z)    ((z).len > 2)
  324. #endif
  325.  
  326.  
  327. /*
  328.  * ztofull - convert a ZVALUE to a FULL if possible
  329.  *
  330.  * If the value is too large, only the low order bits that are able to 
  331.  * be converted into a FULL will be used.
  332.  */
  333. #define ztofull(z)    (zistiny(z) ? ((FULL)((z).v[0])) :        \
  334.                       ((FULL)((z).v[0]) +        \
  335.                        ((FULL)((z).v[1]) << BASEB)))
  336.  
  337. #define z1tol(z)    ((long)((z).v[0]))
  338. #define z2tol(z)    (((long)((z).v[0])) + \
  339.                 (((long)((z).v[1] & MAXHALF)) << BASEB))
  340. #define    zclearval(z)    memset((z).v, 0, (z).len * sizeof(HALF))
  341. #define    zcopyval(z1,z2)    memcpy((z2).v, (z1).v, (z1).len * sizeof(HALF))
  342. #define zquicktrim(z)    {if (((z).len > 1) && ((z).v[(z).len-1] == 0)) \
  343.                 (z).len--;}
  344. #define    zfree(z)    freeh((z).v)
  345.  
  346.  
  347. /*
  348.  * Output modes for numeric displays.
  349.  */
  350. #define MODE_DEFAULT    0
  351. #define MODE_FRAC    1
  352. #define MODE_INT    2
  353. #define MODE_REAL    3
  354. #define MODE_EXP    4
  355. #define MODE_HEX    5
  356. #define MODE_OCTAL    6
  357. #define MODE_BINARY    7
  358. #define MODE_MAX    7
  359.  
  360. #define MODE_INITIAL    MODE_REAL
  361.  
  362.  
  363. /*
  364.  * Output routines for either FILE handles or strings.
  365.  */
  366. extern void math_chr MATH_PROTO((int ch));
  367. extern void math_str MATH_PROTO((char *str));
  368. extern void math_fill MATH_PROTO((char *str, long width));
  369. extern void math_flush MATH_PROTO((void));
  370. extern void math_divertio MATH_PROTO((void));
  371. extern void math_cleardiversions MATH_PROTO((void));
  372. extern void math_setfp MATH_PROTO((FILE *fp));
  373. extern char *math_getdivertedio MATH_PROTO((void));
  374. extern int math_setmode MATH_PROTO((int mode));
  375. extern long math_setdigits MATH_PROTO((long digits));
  376.  
  377.  
  378. #ifdef VARARGS
  379. extern void math_fmt();
  380. #else
  381. extern void math_fmt MATH_PROTO((char *, ...));
  382. #endif
  383.  
  384.  
  385. /*
  386.  * The error routine.
  387.  */
  388. #ifdef VARARGS
  389. extern void math_error();
  390. #else
  391. extern void math_error MATH_PROTO((char *, ...));
  392. #endif
  393.  
  394.  
  395. /*
  396.  * constants used often by the arithmetic routines
  397.  */
  398. extern HALF _zeroval_[], _oneval_[], _twoval_[], _tenval_[];
  399. extern ZVALUE _zero_, _one_, _ten_;
  400.  
  401. extern BOOL _math_abort_;    /* nonzero to abort calculations */
  402. extern ZVALUE _tenpowers_[2 * BASEB];    /* table of 10^2^n */
  403. extern int _outmode_;        /* current output mode */
  404. extern LEN _mul2_;        /* size of number to use multiply algorithm 2 */
  405. extern LEN _sq2_;        /* size of number to use square algorithm 2 */
  406. extern LEN _pow2_;        /* size of modulus to use REDC for powers */
  407. extern LEN _redc2_;        /* size of modulus to use REDC algorithm 2 */
  408. extern HALF *bitmask;        /* bit rotation, norm 0 */
  409.  
  410. #endif
  411.  
  412. /* END CODE */
  413.